r 


AD-769  893 

CALCULATIONS  OF  NEAR  FIELD  PARTHniiAirc 
GROUND  MOTION  EARTHQUAKE 


J.  Theodore  Cherry 
Systems,  Science  and  Software 


[ 


Prepared  for: 

Air  Force  Office  of  Scientific  Research 


29  June  1973 


DISTRIBUTED  BY: 


Nation?!  Technical  Information  Service 
U.  S.  DEPARTMENT  OF  COMMERCE 

5?85  Port  Royal  Road,  Springfield  Va.  22151 


rDOCL\.',i:KY  COMTIiOl.  DATA  -  ft  L  D 

(Security  m  t/i*  si  t  ie  btinn  o  (  f/f/r,  bodv  ct  nbxtrmet  anti  it'.'  n  ir.otution  pxisl  he  entered  when  the  over  *11  report  Is 


UNCLASSIFIED 

Seruntv  (  j i'ifu-nnoti 


l  Ofl  'OINA  TiNE  TIVITV  (Cvrpnrmte  euthor)  Ui.  Rt^CHT  StUr*iTv  CLASSIFICATION 

Systems,  Science  ana  Software  UNCLASSIFIED 

PO  Box  1620  ,t  CDOUP 

La  Jolla,  CA  92037 


3  REPORT  TITLE 

Calculations  of  Near  Field,  Earthquake  Ground  Motion 


4  OESCRiPTlvt  NOTES  (Type  al  report  end  Inclusive  dmte*) 

Scientific . Interim,  Annual.  Technical  Report  1  May  1972 


1*  A  JTMORU  I  (J’tfut  n*«n.  e.  /n/c-je  Int  1 i«  *,  to  st  name) 

J .  T  Cherry 


31  May  1273 


6  RE  FOR  T  O  A  T  C 

29  June  1973 


M.  CONTRACT  ON  CHANT  NO 

F44620-72-C-0051 

b.  PKOJEC  1  NO 

AO  2134 


7m.  TOTAL  NO.  Q  F  PACES  |76.  NO-  OF  REFS 


62701E 


vb.  OTHfR  M*  PORT  NOiSt  (Any  other  numberm  thmt  a.ny  be  ex  si 
thlm  report » 

NrQSR  r  7ft-  73  •  x  TZ 


10  OISTRIBUTION  STATEMENT 


Approved  for  public  release;  distribution  unlimited 


i  «  UjDO'  rwE‘i  »ir*y  N^^ES 


TECH,  OTHER 


ij.  •  *.v  *  -  v:t» 

Air  Force  Office  of  Scientific  Research/'. 
1400  Wi'so.i  Boulevard 
Arlingt  a,  VA  22201 


11-  AOS  TRAC  T 


A  stick-slip  rupture  model  has  been  incorporated  into  a  two-dimensional  (plane 
strain)  Lagrangian  code  with  the  added  feature  that  rupture  initiation  is  plastic 
work  dependent.  Realistic  laboratory  test  data  have  been  used  as  input  to  the 
model.  Theoretical  seismograms  have  been  generated  in  both  the  frequency  and  time 
domain  along  with  a  decomposition  of  the  ground  motion  into  P  and  S  components. 

The  calculations  suggest  how  estimates  of  fault  length,  rupture  velocity  and 
dynamic  stress  drop  may  be  obtained  from  near  field  data.  These  estimates  provide 
an  initial  specification  of  the  parameters  required  by  the  rupture  model  in  order 
to  match  a  given  set  of  near  field  data. 


Reproduced  by 

NATIONAL  TECHNICAL 
INFORMATION  SERVICE 

U  S  Deportment  of  Con  nerct 
Springfield  VA  27.!51 


UNCLASSIFIED 


AD  76  9893 


•HOSH-TR- 


78  -  lia1* 


SYSTEMS,  SCIENCE  AND  SOFTWARE 


SSS-R-75-1759 


CALCULATION  OF  NEAR  FIELD,  EARTHQUAKE  GROUND  MOTION 


Annual  Technical  Report 
For  Period  Ending  31  May  1973 


Project  Manager:  J.  Theodore  Cherry  (714)  453-0060 


Sponsored  by 

Advanced  Research  Projects  Agency 
ARPA  Order  No.'  2134 


Monitored  by 

Air  Force  Office  of  Scientific  Research 

Arlington,  Virginia  22209  [[) 

fnV?^ 


n 


Contract  No.  F44620-72-C-0051 
Program  Code  2F10 

Effective  Date  of  Contract:  1  June  1972 
Contract  Expiration  Date:  28  Feb.  1974 
Amount  of  Contract:  $112,709 


C 

1  FRF' 


NOV  14  1973 

IbiiSKO  U  IS 

B 


S  Project  195 


j  Kfmi j&Htsr 

EritelbnVm 


June  29,  1973 


P  O.  BOX  1620,  LA  JOLLA,  CALIFORNIA  92037,  TELEPHONE  (714)  453-0060 


) 


» 

SYSTEMS,  SCIENCE  AND  SOFTWARE 


> 


sss-R-73-r/sy 


> 


> 


CALCULATIONS  OF  NEAR  FIELD,  EARTHQUAKE  GROUND  MOTION 

Annual  Technical  Report 
For  Period  Ending  31  May  1973 

Project  Manager:  J.  Theodore  Cherry  (714)  453-0060 


Sponsored  by 

Advanced  Research  Projects  Agency 
ARPA  Order  No.  2134 


Monitored  by 

Air  Force  Office  of  Scientific  Research 
Arlington,  Virginia  22209 


Contract  No.  F4462 
Program  Code  2F10 
Effective  Date  of 
Contract  Expiratio 
Amount  of  Contract 


June  29,  1973 


^-72-C-0051 

Contract:  1  June  1972 
n  Date:  28  Feb.  1974 
:  $112,709 


r 


Dio.  - 
Appr  -  ‘  • 


juA 

,jC  rMcasoi 

-vU*d 


S3  Project  195 


P.  O.  BOX  1620,  LA  JOLLA,  CALIFORNIA  92037,  TELEPHONE  (714)  453-0060 


4 


FOREWORD 


This  annual  technical  report  entitled  "Calculations 
of  Near  Field,  Earthquake  Ground  Motion,"  is  submitted  by 
Systems,  Science  and  Software  (S3)  to  the  Advanced  Research 
Projects  Agency  (ARPA)  and  to  the  Air  Force  Office  of 
Scientific  Research  (AFOSR) . 

The  report  presents  the  results  obtained  over  the 
first  twelve  months  of  an  ei ghteen-month  effort  designed  to 
uncover  the  relation  between  earthquake  ground  motion  and 
the  nature  of  the  earthquake  source. 

This  research  was  supported  by  the  Advance..  Research 
Projects  Agency  of  the  Department  of  Defense  and  was  monitored 
by  the  Air  Force  Office  of  Scientific  Research  under  Contract 
No.  F44620-72-C-0051 .  Mr.  Donald  C.  Clements  was  the  ARPA 
Program  Manager  and  Col.  Donald  V.\  IClick  the  AFOSR  Project 
Scientist. 

Dr.  J.  T.  Cherry  was  the  S3  Project  Manager  for  the 
study.  Mr.  E.  J.  Ilalda  was  responsible  for  the  scientific 
programming  associated  with  the  ground  motion  calculations. 

Mr.  K.  G.  Hamilton  developed  fne  computer  codes  required 
to  process  and  display  the  results. 


ABSTRACT 


A  stick-slip  rupture  model  has  been  incorporated  into  a 
two-dimensional  (plane  strain)  Lagrangian  code  with  the  added 
feature  that  rupture  initiation  is  plastic  work  dependent. 

Realistic  laboratory  test  data  have  been  used  as 
input  to  the  model.  Theoretical  seismograms  have  been  gene¬ 
rated  in  both  the  frequency  and  time  domain  along  with  a 
decomposition  of  the  ground  motion  into  P  and  S  components. 

The  calculations  suggest  how  estimates  of  fault  length, 
rupture  velocity  and  dynamic  stress  drop  may  be  obtained  from 
near  field  data.  These  estimates  provide  an  initial  speci¬ 
fication  of  the  parameters  required  by  the  rupture  model  in 
order  to  match  a  given  set  of  near  field  data. 
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I.  INTRODUCTION 

The  elastic  rebound  theory,  developed  by  H.  F.  Reid^ 
after  the  1906  San  Francisco  earthquake  identifies  the  imme¬ 
diate  source  of  an  earthquake  as  the  release  of  accumulated 
strain  energy  in  the  rock  mass  surrounding  a  fault.  This 
statement  implies  that  the  rupture  process  at  the  fault  sur¬ 
face  is  responsible  for  the  relative  displacement  across  the 
fault  and  hence  the  seismic  waves  radiated  during  the  eartnquake. 
Almost  sixty  years  passed  before  Benioff1  J  showed  that  the 
displacement  at  the  fault  inferred  from  the  elastic  rebound 
theory  was  consistent  with  seismological  observations  which 
suggested  that  the  earthquake  source  function,  for  the  radiated 
seismic  energy,  should  be  a  double  couple. 

Following  Benioff's  work  it  became  clear  that  an  under¬ 
standing  of  the  nature  of  the  earthquake  source  depended  on 
an  adequate  simulation  of  fracture  propagation  over  the  fault 
plane.  Analytic  techniques  were  applied  to  the  problem. 

Analytic  models  have  assumed  either  a  moving  disloca¬ 
tion^3,  ^  or  a  prescribed  stress  relaxation ^ order  to 
simulate  the  propagating  rupture.  Differences  between  these 
techniques  should  occur  only  because  the  assumed  stress  relaxa¬ 
tion  does  not  produce  the  same  time  history  of  relative  displace¬ 
ment  at  the  fault  surface  as  that  used  during  a  specific  exercise 
of  the  dislocation  technique. 

Obviously  these  analytic  models  assume  the  behavior  of 
the  equivalent  elastic  source  at  the  fault  surface.  This 
assumption  is  stated  either  in  a  dislocation  or  stress  relaxa¬ 
tion  format.  If  analytic  techniques  are  to  be  used  with  confi¬ 
dence  then  arbitrary  assumptions,  regarding  the  behavior  of 
the  equivalent  elastic  source  at  the  fault  surface,  must  be 
removed . 
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When  analytic  models  are  used  to  match  the  free  field 
ground  motion  from  a  specific  earthquake  in  the  region  where 
the  material  response  is  linear-elastic,  only  then  will  the 
equivalent  elastic  source  for  the  earthquake  become  available. 

The  match  implies  tnat  the  analytic  technique  has  used  the 
correct  description  of  the  equivalent  source  at  the  fault 
surface . * 

Free-f ield  ground  motion  measurement  from  actual  earth¬ 
quakes  do  not  exist.  A  few  free  surface,  iu..r  field  measurements 
have  been  made  (for  the  Parkfield,  June  28,  1966,  and  San  Fernando, 
February  9,  1971,  earthquakes)  and  more  data  of  this  type  are 
certain  to  become  available.  It  is  important,  therefore,  that 
techniques  be  developed  which  are  capable  of: 

1.  Simulating  the  rupture  process  in  such  a  manner 
that  laboratory  data  from  appropriate  rock  tests 
may  be  used  to  specify  the  parameters  in  the 
rupture  model. 

2.  calculating  the  theoretical,  free-field  seismograms 
caused  by  the  rupture  in  order  to  obtain  the  equi¬ 
valent  elastic  source. 

3.  Including  the  effect  of  the  free  surface  and  local 
site  geology  in  order  to  compare  calculated  ground 
motion  with  free- surface ,  close-in  measurements. 

Lagrangian  computer  codes  have  been  developed ^ 
which  are  capable  of  simulating  the  response  of  geologic 
materials  to  a  propagating  stress  wave  of  arbitrary  amplitude. 

*rrr°urse  acbual  behavior  of  the  fault  surface  will  be  quite 
different  from  tnat  specified  by  the  equivalent  source  if  non¬ 
linear  material  behavior  occurs  during  the  rupture  process. 

An  analogous  situation  is  encountered  for  spherical  explosions 
j-rr  ac^ual  material  response  in  the  nonlinear  region 

is  different  from  that  obtained  using  the  equivalent  elastic 
source  (the  reduced  displacement  potential). 
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These  codes  are  capable  of  carrying  the  calculations  into  the 
small  displacement,  elastic  region  and  yet  flexible  enough 
to  permit  very  general  material  response  formulations  in  the 
nonlinear  region.  ^ ^  They  have  been  used  extensively  to 
both  predict  the  effects  of  explosive  sources  on  the  surrounding 
rock  environment ^ and  to  obtain  the  equivalent  elastic 
source  as  a  function  of  rock  type,  depth  of  burial,  and  explo¬ 
sive  yield. 

? 

In  an  attempt  to  at  least  partially  satisfy  the  above 
three  objectives,  a  stick-slip  rupture  model  has  been  incorporated 
into  a  two-dimensional  (plane  strain)  Lagrangian  stress  wave 
f  code.  This  earthquake  model  now  furnishes  the  near  source 

(free-field  and  free  surface)  ground  motion  caused  by  the  stick- 
slip  rupture  process.  The  only  limitation  in  the  model  is  the 
plane  strain  assumption.  This  assumption  implies  an  infinite 
t  fault  dimension  normal  to  the  plane.  If  the  calculations  were 

carried  to  their  far-field  limit,  they  would  correspond  to  a 
line  source  rather  than  a  point  source.* 

In  spite  of  this  assumotion,  this  two-dimensional  fault 
7 

model  has  provided  a  great  deal  of  insight  into  the  near  field 
variation  of  peak  particle  velocity,  corner  frequency  and 
displacement  spectra  with  fault  length,  rupture  velocity  and 
stress  drop.  The  physics  of  the  source  seems  to  be  easy  to 
simulate  (certainly  easier  than  nuclear  explosions).  The 
last  remaining  obstacle  to  be  overcome,  before  teleseismic 
amplitude  dependence  can  be  related  to  the  nature  of  the 
earthquake  source,  is  the  inclusion  of  the  third  space  dimen¬ 
sion  in  the  model. 


Theoretical  seismograms  from  this  model  would  be  appropriate 
for  an  actual  earthquake  having  an  out-of-plane  fault  dimen- 
I  sion  comparable  to  the  distance  between  the  fault  and  the 

seismome  ter . 
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II.  the  physics  of  the  earthqu/vke  source 


A  formulation  of  the  rupture  process  must  provide  quan¬ 
titative  answers  to  the  following  three  questions: 

1.  Why  does  rupture  occur? 

2.  What  is  the  stress  adjustment  during  rupture? 

j.  When  does  the  rupture  heal? 

A  stick- slip  model  of  rupture  has  been  formulated  which  answers 
the  above  questions  as  follows: 

1.  Rupture  initiation  is  plastic  work  dependent. 

2.  During  rupture  the  tangential  stress  at  the 
slipping  interface  is  relaxed  to  its  kinetic 
friction  value.  This  relaxation  allows 
adjacent  points  on  the  interface  to  move 
apart  (slip) . 

3.  The  rupture  heals  (adjacent  points  on  the  inter¬ 
face  stick)  if  the  relative  velocity  of  two 
adjacent  points  changes  sign  and  if  the  tangential 
stresses  at  the  interface  is  sufficient  to  main¬ 
tain  continuity  of  tangential  velocity. 

A  Lagrangian  code  is  ideally  suited  to  simulate  a 
slipping  interface.  It  is  simply  necessary  to  first  decouple 
the  grid  line,  over  which  slip  is  to  occur,  in  order  to  isolate 
the  normal  and  tangential  components  of  stress  at  the  inter¬ 
face  (Equations  22  through  25,  Appendix  I)  and  second,  to 
apply  contact  discontinuity  boundary  conditions,  involving 
continuity  of  normal  stress  and  normal  velocity  components, 
in  order  to  solve  for  the  normal  stress  component  (Eq.  18, 
Appendix  II).  If  the  boundary  point  is  "welded"  (not  slipping) 
then  the  tangential  velocity  component  will  also  be  continuous. 
This  latter  condition  permits  a  unique  solution  for  the 
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tangential  stress  at  the  welded  point  (Eq.  19,  Appendix  II). 

The  basic  mechanism  for  releasing  the  strain  energy  in  this 
rupture  model,  is  the  relaxation  of  the  tangential  stress  from 
its  "welded"  value  to  its  "kinetic  friction"  value. 

Figure  1  follows  the  tangential  stress  at  two  points  on 
the  fault  surface  (slipping  interface)  during  a  free-field 
calculation  in  which  the  stick-slip  model  was  used  to  release 
the  strain  energy.  For  this  particular  calculation  the  fault 
length  evf  ntually  grew  to  10  km.  The  solid  curve  in  the  figure 
corresponds  to  the  point  on  the  fault  where  the  rupture  starts 
(the  focus)  while  the  dashed  curve  is  for  a  point  on  the  fault 
2.5  km  away.  The  initial  value  of  tangential  stress  (t  )  on 
the  fault  was  o xe  kbar.  During  rupture  this  stress  component 
is  relaxed  to  its  kinetic  friction  value  (x^)  .  In  this  problem 
rk  =  0.5  kbar.  The  tangential  stress  is  maintained  at  the  rk 
value  until  adjacent  points  on  each  side  of  the  fault  reverse 
velocity.  When  the  velocity  reversal  occurs,  the  points  are 
tied  (the  fault  sticks)  if  the  tangential  stress,  required  to 
maintain  continuity  of  tangential  velocity  lies  between  +  rk. 
After  the  points  are  tied  the  tangential  stress  finds  a  static 
equilibrium  value  (rs). 

Notice  that  for  the  point  2.5  km  away  from  the  beginning 
of  the  fault,  the  tangential  stress  builds  to  a  maximum  of 
1.43  kbar  due  to  stress  differences  parallel  to  the  fault,  be¬ 
fore  the  plastic  work  criterion  at  this  distance  is  violated. 
This  occurs  a;  1.12  seconds;  rupture  begins  and  the  tangential 
stress  is  relaxed  to  Tk. 

1  igure  2  shows  the  final  static  level  attained  by  the 
tangential  stress  over  the  10-km  fault.  Results  from  all 
problems  run  to  date  indicate  that  at  least  one  point  on  the 
fault  will  stick  early  and  cause  a  mild  stress  concentration 
to  occur  in  the  static  solution.  In  this  figure  the  concentra¬ 
tion  occurs  3  km  from  the  focus.  Figure  2  also  shows  that  most 
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Fig.  1- -Tangential  stress  versus  time  at  two  points  2.5  km  apart  on  the  fault 
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of  the  static  stress  drop  occurs  at  the  end  of  the  fault  where 
fhe  rupture  stops;  a  result  that  is  again  common  to  all  cal¬ 
culations  having  a  finite  rupture  velocity. 

*  Figure  3  shows  the  relative  displacement  over  the  fault 
at  0. '-second  intervals.  The  point  that  ties  early  at  3  km  is 
responsible  for  the  stress  concentration  at  that  distance  in 
Fig.  2. 

f 

The  rupture  velocity  over  the  first  7  km  of  the  fault 
was  2.15  km/sec.  At  this  distance  the  plastic  work  criterion 
was  increased  so  that  the  rupture  would  not  stop  abruptly, 
f  Figure  4  shows  tin  arrival  time  of  the  rupture  versus  distance 

along  the  fault.  Over  the  last-  3  km  the  rupture  velocity  is 
approximately  1.6  km/sec,  giving  an  average  rupture  velocity 
over  the  entire  fault  of  2  km/sec. 

*  Plastic  flow  is  due  to  the  inability  of  real  geologic 
materials  to  support  unlimited  values  of  shear  stress.  The 
deviatoric  stress  components  in  the  yielding  element  are 
modified  such  that  the  resulting  stress  state  is  consistent 
with  a  Mises  yield  criterion  (Eq.  1,  Appendix  III). 

Initially,  plastic  flow  was  included  in  the  model  in 
order  to  remove  the  large  stress  concentrations  that  occurred 
I  the  ends  of  the  rupture  during  calculations  involving  only 

linear,  elastic  material  behavior.  It  was  immediately  found 
that  rupture  velocity  could  be  controlled  by  allowing  rupture 
initiation  to  be  dependent  on  the  plastic  work  dissipated 
I  during  the  yielding  process  (Eqs .  3  th-ough  5,  Appendix  III). 

There  is  experimental  evidence that  crystalline 
rocks  undergo  significant  yielding  prior  to  brittle  failure 
at  the  temperatures  and  pressures  appropriate  even  for  shallow 
t  earthquakes  (focal  depths  of  around  10  km).  We  have  allowed 

this  mechanism  to  control  the  rupture  velocity  by  specifying 
the  plastic  work  for  rupture  to  be  a  function  of  distance  from 
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the  focus  (Eq.  6,  Appendix  III).  This  requires  not  only  that 
the  dimensions  of  the  fault  zone  be  specified,  but  also  the 
relation  between  the  yield  surface  and  the  stress  state  in 
the  fault  zone. 

The  calculations  reported  in  the  next  section  were  ob¬ 
tained  assuming  an  elliptical  fault  zone  with  all  the  material 
in  the  fault  zone  initially  lying  on  the  yield  surface,  i.e., 
an  attempt  to  increase  the  second  deviatoric  invariant  above 
its  initial  value  in  the  fault  zone  causes  plastic  flow  and 
therefore  plastic  work. 

The  code,  of  course,  is  not  limited  to  a  plastic  work 
rupture  criterion.  In  the  code,  the  components  of  stress  have 
been  isolated  at  the  fault  surface  and  a  rupture  criterion 
could  easily  be  formulated  in  terms  of  these  stress  components. 

For  example,  the  tangential  stress  could  initially  be 
limited  at  the  boundary.  This  would  allow  the  fault  to  slide 
stabl> ,  i.e.,  creep.  (In  Fig.  1,  a  creep  event  would  occur 
at  2.5  km  of  the  allowable  tangential  stress  at  this  distance 
was  less  than  1.43  kbar).  The  drop  in  tangential  stress  to 
its  kinetic  friction  value  (rupture)  could  then  be  made  a 
function  of  the  size  of  the  creep  event.  Stable  sliding  has 
been  observed,  prior  to  rupture,  in  laboratory  stick-slip 
events  in  plates  of  Westerly  granite. Creep  may  prepare 
the  fault  surface  for  rupture  by  polishing  the  surface. 

Rupture  velocity  could  be  controlled  by  varying  the  magnitude 
of  the  creep  event  required  to  cause  the  tangential  stress 
to  drop  to  its  kinetic  friction  value. 

Frictional  sliding  on  ground  surfaces  of  granite  has 
been  investigated  by  Byerlee.^9^  For  normal  stresses  (a  ) 
varying  between  2  -  12  kbar  over  the  surface,  he  found  that 
the  tangential  stress  drop  from  static  friction  (t  )  to 
kinetic  friction  (xk)  is  given  by 
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To  ’  Tk  =  0,25  +  0,13  an  (a11  units  in  kbars)  (2.1) 

f  While  this  relation  was  established  for  specimens  at  room 

temperature  and  fo’~  specially  prepared  surfaces,  it  probably 
furnishes  an  upper  limit  to  the  allowable  d)namic  stress 
drop  for  shallow  earthquakes.  Laboratory  experiments  should 
1  performed  in  order  to  (1)  determine  che  effect  of  tempera- 

tue  on  dynamic  stress  drop  and  (2)  to  determine  the  amount  of 
inelastic  work  required  to  cause  stick-slip.  Results  of  these 
experiments  would  provide  the  input  required  by  the  stick-slip 
.  rupture  model. 
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III.  CALCULATED  GROUND  MOTION  RESULTING  FROM  THE 
STICK-SLIP  RUPTURE  MODEL 

3.1  SUMMARY  OF  INPUT  PARAMETERS 


In  order  to  exercise  the  rupture  model  in  its  current 
form,  a  number  of  parameters  must  be  specified.  These  are: 

9  The  Fault  Zone.  An  elliptical  fault  zone,  defined 
as 


r 


t 


t 


f 


» 


(3.1) 


was  assumed  with  all  the  material  in  the  fault  zone 
initially  lying  on  the  yield  surface.  The  origin 
of  the  x-y  coordinate  system  is  located  at  the  center 
of  the  fault,  as  shown  in  Fig.  5.  Elastic  behavior 
was  assumed  for  the  material  outside  the  fault  zone. 
The  minor  axis,  b,  was  2  km  for  all  calculations. 

The  major  axis,  a,  extended  3  km  beyond  the  end 
of  the  fault. 


•  The  Plastic  Work  Required  for  Pupture  Initiation. 

A  unidirectional  rupture  was  assumed  as  shown  in 
Fig.  5.  For  a  given  fault  length,  L,  the  functional 
form  used  to  initiate  rupture  within  the  interval 
-L/2  <_  x  <_L/2  was  assumed  to  be  (Appendix  III, 

Eq.  6) 


W 


1  1  x+L/2 

7  '  T 


(°  £  x  +  I  1  f)  (3.2a) 


1  7  x+L/2 

‘  7  3  ~T~_ 


(x  *  T  -  l)  (3'2b) 
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Rupture  starts  at  x  =  -L/2  (W  =  0)  and  is  assumed  to  terminate 
at  x  =  L/2.  In  the  calculations,  rupturing  was  not  permitted 
outside  the  interval  -L/2  ^  x  £  L/2. 

•  The  Stress  Drop  During  Rupture  (t  -  x,J  .  It  was 

.  o  K 

found  that  the  normal  stress  (a  )  on  the  fault 

n 

did  not  change  from  its  initial  value  during  the 
calculation.  Therefore,  the  assumed  stress  drop 
remained  constant  over  the  fault.  Equation  (2.1) 
was  used  as  a  guide  for  this  quantity  and  calcu¬ 
lations  were  obtained  for  stress  drops  of  0.25  and 

0.5  kbars.  The  initial  stress  (t  )  was  the  ca;... 

o 

for  all  calculations  with  t  =1  kbar. 

o 

•  The  Shes r  Modulus  (p) ,  Bulk  Modulus  (k)  and  Density  (p) . 
These  quantities  were  maintained  constant  for  all 
calculations,  with  p  =  324  kbar,  k  =  478  kbar  and 

p  =  2.8  g/cc.  The  corresponding  compression  (a) 
and  shear  (8)  wave  velocities  are  a  =  5.7  km/'sec, 

8  =  3.4  km/sec. 

Table  I  summarizes  the  input  parameters  (a,  c,  d,  L, 
t  -  T-i  )  along  with  the  rupture  velocity  (V  and  static  stress 

0  K 

drop  (t  -  t  )  that  resulted  from  the  four  calculations.  Since 
o 

a,  8,  p,  b,  and  t  were  the  same  for  all  calculations,  these 

o 

parameters  are  not  listed  in  the  table. 


t 


> 


R- 1759 


I 


t 


f 


TABLE  I 


Calculation 

a  (km) 

c 

d  (km) 

L(km) 

T  -T, 
o  k 

(kbar) 

VR 

(km/sec) 

t  -t 

0  S 

(kbar) 

10A 

8 

0.7 

10 

10 

0.5 

2.0 

0.146 

5A 

5.5 

0.7 

10 

5 

0.5 

2.15 

0.178 

5B 

5.5 

0.4 

10 

5 

0.5 

3.75 

0.224 

5C 

5.5 

0.4 

10 

5 

0.25 

2.15 

0.110 

1 
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3 . 2  RADIATION  PATTERNS 

All  ground  motion  calculations  are  easily  separable 
into  P  and  S  components.  This  is  accomplished  by  monitor- 

-y. 

ing  the  divergence  (V  •  s)  and  curl  (V  x  s)  of  the  dis¬ 
placement  field  at  selected  points  in  the  elastic  regime. 

The  equation  of  motion  for  a  linear  elastic  stress 

wave  is 


r 


t 


t 


r 


f 


f 


» 


—  =  a2  7  (V  •  s)-32  Vx(V  x  s)  (3.3) 

3t 2 

where  s  is  particle  displacement  and  a  and  3  are  the 
compression  and  shear  velocities.  If  s  is  separated  into 
a  scalar  (4>)  and  vector  (x)  potential,  then 

s  =  -V<f.  -  Vxy  (V  •  x  =  0)  (3.4) 

Substituting  Eq.  (3.4)  into  Eq.  (3.3)  gives  the  scalar  and 
vector  wave  equations 

^  =  a2  V2<J>  (3.5) 

3 1 2 

^  =  32  V2x  (3.6) 

3 1 2 

Hooke's  law  gives 

P  =  -k  V  •  s  (3.7) 

where  k  is  the  bulk  modulus,  and  p  is  the  pressure  com¬ 
ponent  of  the  stress  tensor.  Also,  since  V  •  (V  x  y)  =  0, 
then  from  Eq.  (2) 

V  •  s  =  -V2*  (3.8) 
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Using  Eqs .  (3.5),  (3.7)  and  (3.8)  gives 


3  2  4 
3t 2 


a 

TT 


P  ■ 


(3.9) 


where  y  is  th<-  siiear  modulus  and  p  the  density.  Therefore 
<{)  is  easily  obtained  from  a  calcuTation  since  p  is  a 
saved  variable  in  the  code  . 

Sir.^e  v  x  ( V4, )  =  0  and  v  x  (y  x  ^)  =  -V2x 


then  Eq.  (3.4)  gives 

V  x  s  =  V2x  (3.1  ) 

Equations  (3.6)  and  (3.10)  give 

*  2 

— =  62  (V  x  s)  where  B2  =  —  (3.11) 

3 1 2  9 

The  quantity  V  x  s  is  proportional  to  the  rotation  of 
a  line.element  and  is  also  a  saved  variable  in  the  code.  There¬ 
fore  x  is  readily  nvi.lnble  from  the  calculations.  The  form 
of  Eq.  (3.11)  used  in  the  code  is 


x(t) 


dt 


(3.12) 


•  • 

where  x  and  y  are  particle  velocities  in  the  x  and  y 
directions . 

Equations  (3.9)  and  (3.12)  permit  the  calculated  ground 
motion  to  be  separated  into  P  and  S  components.  A  posi¬ 
tive  value  of  <J)  corresponds  to  a  compression  (p  >  0)  while 
a  positive  value  of  x  gives  a  clockwise  body  rotation. 
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Figure  5  shows  the  stations  monitored  during  each 
calculation.  The  origin  of  the  x-y  coordinate  system  is 
located  at  the  center  of  the  fault.  The  rupture  velocity 
was  unidirectional  and  propagated  in  the  positive  x  direc¬ 
tion  . 

Figures  6  through  13  show  the  P($)  and  S(x)  radia¬ 
tion  patterns  from  the  four  calculations.  These  figures  give 
the  peak  values  of  the  two  potentials  at  all  stations  located 
10  krr  r->om  the  center  of  the  fault.  The  finite  rupture 
velocity  produces  a  noticable  Doppler  effect  in  the  radia¬ 
tion  patterns,  with  the  S  component  n  the  fault  plane 
(in  the  rupture  direction)  approximately  2.5  times  larger 
than  the  S  components  in  the  auxiliary  plane.  The  same 
factor  applies  to  the  smar/pmax  ratio. 

While  the  value  for  Snax  in  the  auxiliary  plane 
is  real,  the  peaks  in  the  other  lobes  in  both  P  and  S 
may  have  been  missed.  This  is  due  to  an  inadequate  azi¬ 
muthal  sampling  along  the  10-km  radius  (Fig.  5).  While 
this  will  be  remedied  in  future  calculations,  it  appears 
that  the  routine  for  P  and  S  separation  produces 
acceptable  radiation  patterns. 

Figures  14  through  21  give  the  full  time  history  of 
<j)  and  x  at  Station  2  for  the  four  calculations  and  show 
the  basic  difference  between  an  earthquake  source  and  a 
center  of  dilatation  (an  explosion).  The  center  of 
dilatation  gives  a  scalar  potential  (<{>)  that  resembles 
4>  for  an  earthquake.  Of  course  the  vector  potential 
vanishes  for  the  center  of  dilatation.  The  calculations 
show  that  the  earthquake  source  produces  static,  aximuthally 
dependent  values  of  pressure  (Eq.  (3.9))  and  body  rotation 
(Em  (3.12)).  This  feature  should  cause  an  increase  in 
surface  wave  excitation  for  earthquakes  relative  to  explosions 
and  may  be  responsible  for  the  effectiveness  of  the  Ms/mb 
discriminant . 
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Fig.  9--S  radiation 
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Fig.  15--Time  history  of  vector  potential  (x)  at  Station 
2,  10A. 
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Fig.  17--Time  history  of  vector  potential  (x)  at  Station 
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Fig.  20--Timc  history  of  scalar  potential  (4>)  at  Station 
2,  5C . 
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3.3  PEAK  PARTICLE  VELOCITY 

Peak  particle  velocity  was  determined  at  all  stations 
10  km  from  the  center  of  the  fault.  Figure  22  shows  peak 
particle  velocity  plotted  versus  the  product  of  rupture 
velocity  <VR)  and  dynamic  stress  drop  (t  -  ij,)  .  The 
functional  relation  seems  to  be  linear,  with  the  slope 
varying  with  azimuth  as  shown.  Fault  length  does  not  seem 
to  be  important  in  determining  peak  velocity. 

Figures  23  through  42  show  the  velocity  seismograms 
for  the  four  calculations  at  Stations  1  through  5.  The 
complexity  of  the  seismogram  increases  with  increasing  azimuth, 
apparently  due  to  the  temporal  separation  of  the  starting  and 
stopping  phases  of  the  propagating  rupture. 
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Fig.  23--Particle  velocity  at  Station  1,  10A. 
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Fig.  24  -  -  Particle  velocity  at  Station  2,  10A. 
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Fig.  25--Particle  velocity  at  Station  3,  10A. 
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Fig.  27--Particle  velocity  at  Station  5,  10A. 
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Fig.  28--Particle  velocity  at  Station  1,  5A. 
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Fig.  33- -Particle  velocity  at  Station  1,  5B. 
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Fig.  35--Particle  velocity  at  Station  3,  5B. 
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Fig.  37--Particle  velocity  at  Station  5,  5B. 
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3.4  DISPLACEMENT  SPECTRA  AND  CORNER  FREQUENCY 

i  ;  7 - * - 

A  sketch  of  a  typical  near  field  displacement  spectrum 
is  shown  in  Fig.  43.  At  low  frequencies  the  spectrum  shows 
an  f- 1  trend,  caused  by  the  static  displacement.  At  high 

t  frequencies  the  spectrum  decays  approximately  as  f * 3 ,  in¬ 

dicating  that  the  equivalent  elastic  source  should  be  con¬ 
tinuous  in  both  displacement  and  particle  velocity. 

As  shown  in  Fig.  43,  the  corner  frequency  is  defined 
to  occur  at  the  transition  to  the  high  frequency  (f*3)  trend. 
Figure  44  shows  corner  frequency  plotted  versus  VD/L  for  the 
indicated  stations  at  a  radius  of  10  km  from  the  center  of 
the  fault.  Corner  frequency  seems  to  be  linearly  related  to 
the  Vp/L  ratio,  with  the  slope  varying  with  a; imuth  as  shown. 
Apparently,  corner  frequency  is  independent  of  dynamic 
stress  drop. 

9  Figure  45  through  48  show  the  x  and  y  components 

of  displacement  spectra  at  Station  2.  The  corner  frequencies 
are  marked  in  each  figure.  A  less  ambiguous  determination  of 
the  corner  frequency  should  be  possible  when  the  near  and  far 
field  source  components  are  separated. 

*«  * 

Finally,  there  is  some  indication  from  the  spectra  shown 
in  Figs.  45  through  48  that  the  value  of  the  spectrum  at  the 
corner  frequency  is  proportional  to  the  product  of  the  fault 

1  length  (L)  and  the  dynamic  stress  drop  (x  -  x,  ) . 
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Fig.  4 
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.  fe 

Frequency 

3- -Sketch  of  typical  near  field  displacement 
spectrum  showing  f" 1  trend  at  low  fre¬ 
quencies,  f " 3  trend  at  high  frequencies 
and  definition  of  corner  frequency  (fc). 
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Fig.  44--Corner  frequency  versus  the  ratio  of  rupture 
velocity  and  fault  length  at  a  radius  of 
10  km  from  the  center  of  the  fault. 
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Fig.  47 - -Displacement  spectra  at  Station  2,  5B. 
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IV.  SUMMARY 
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A  stick- slip  rupture  model  has  been  incorporated  into 
a  two-dimensional  Lagrangian  code  with  the  added  feature  that 
rupture  initiation  is  plastic  work  dependent.  Laboratory 
data,  from  appropriate  rock  mechanics  tests,  may  be  used  to 
specify  the  parameters  in  the  model. 

Theoretical  seismograms  have  been  generated  in  both  the 
frequency  and  time  domain  along  with  a  decomposition  of  the 
ground  motion  into  P  and  S  components.  Results  from  the 
four  calculations,  run  to  date,  are  as  follows: 


A  finite  rupture  velocity  produces  a  noticable 
Doppler  effect  in  the  radiation  pattern  with 
the  S  component  in  the  fault  plane  approxi¬ 
mately  2.5  times  larger  than  the  S  component 
in  the  auxiliary  plane.  The  same  factor  applies 

to  the  S  /P  ratio, 
max  max 

The  earthquake  source  produces  static,  azimuthally 
dependent  values  of  pressure  and  body  rotation. 

The  feature  should  cause  an  increase  in  surface 
wave  excitation  for  earthquakes  relative  to  ex¬ 
plosions  and  may  be  responsible  for  the  effective¬ 
ness  of  the  ^s/mb  discriminant. 

Peak  particle  velocity  (v)  is  linearly  related  to 
the  product  of  rupture  velocity  (V^)  and  dynamic 


stress  drop  (x 


-  V 


From  Fig.  22 


v  =  0.35  Vr(To  -  Tk) 


(4.1) 


where  the  constant  has  units  of  (m) (km) " 1 (kb) ’ 1 
and  represents  an  azimuthal  average. 
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•  At  low  frequencies,  the  displacement  spectrum  shows 
an  f’1  trend,  caused  by  the  local  static  displace¬ 
ment.  At  high  frequencies  the  spectrum  decays  as 
f-3,  indicating  that  the  equivalent  elastic  source 
should  be  continuous  in  both  displacement  and  particle 
velocity . 

•  Corner  frequency  (f  )  seems  to  be  linearly  related 
to  the  ratio  of  rupture  velocity  (V^)  and  fault 
length  (L) .  From  Fig.  44 


f 


c 


(4.2) 


•  There  is  some  indication  that  value  of  the  displace- 

A 

ment  spectrum  (u)  at  the  corner  frequency  is  pro¬ 
portional  to  the  product  of  fault  length  (L)  and 
stress  drop  (t  -  x.) .  From  Figs.  15  through  48 

0  K 

u  =  0 . 4  L  (t  -  t -i  )  (4.3) 

0 

where  the  coefficient  has  units  of  (m) (sec) (km) " 1 (kb) " 1 . 

Equations  (4.1),  (4.2),  and  (4.3)  do  not  permit  a  unique 
determination  of  V^,  L  and  x^  -  x^  given  v,  u  and  fc  in 
the  near  field.  However,  if  limits  can  be  placed  on  one  of  the 
parameters,  the  stress  drop  for  instance,  then  the  remaining 
parameters  may  also  be  limited. 

For  example,  suppose  the  following  near  field  measure¬ 
ments  have  been  obtained: 


f  =  1  HZ 
c 

A 

u  =  8  x  10'2  m-sec  (at  10  km) 
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and  the  stress  drop  (x  -  TjJ  is  assumed  to  be  0.2  kbar. 
Equation  (4.3)  gives 


L  =  1  km 


and  Eq.  (4.2)  gives 


=  2.5  km/sec 

This  implies  that  a  peak  particle  velocity  of  approximately 
18  cm/sec  has  also  been  measured.  If  the  stress  drop  is 
doubled  then  both  the  rupture  velocity  and  fault  length  must 
be  decreased  by  the  same  factor. 

Equations  (4.1),  (4.2)  and  (4.3)  should  be  regarded  as 
a  means  of  obtaining  at  least  some  insight  into  the  parameters 
required  by  the  rupture  model  in  order  to  match  a  given  set  of 
near  field  data.  Iteration  on  the  parameters  will  obviously 
be  required  in  order  to  produce  the  optimum  match. 

Finally,  the  code  is  practically  unlimited  in  terms  of 
its  ability  to  accept  a  given  rupture  model.  This  flexibility 
has  been  obtained  by  isolating  the  normal  and  tangential 
stress  components  at  the  fault  surface  and  solving  for  these 
stress  components  using  contact  discontinuity  boundary  condi¬ 
tions.  A  description  of  the  rupture  process  now  depends  only 
on  an  imaginative  interpretation  and  synthesis  of  appropriate 
laboratory  test  data. 
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APPENDIX  I 


The  difference  equations  used  in  CRAM^10^  to  move  an 
interior  point  are  written  such  that  the  boundary  conditions 
for  an  exterior  point  are  obscured.  Since  we  would  like  to 
isolate  the  normal  and  tangential  r tresses  at  an  interior 
interface,  this  requires  that  the  interface  be  treated  as  an 
exterior  line  over  which  the  correct  boundary  stresses  are 
applied.  In  this  section  a  differencing  scheme  is  obtained 
that  isolates  the  boundary  stresses  and  that  is  consistent 
with  the  CRAM  interior  difference  equations. 

7  The  conservation  of  linear  momentum  (equation  of  motion) 

in  two-dimensional  Cartesian  geometry  is 


CD 

(2) 


If  a  Lagrangian  coordinate  system  (k,j)  is  established  in 
the  material,  then 


3 x  _  ax  ax  ,  ax  ay 
Tic  ax  Tfc  7y  ak 


ax  _  ax  ax  ,  ax  ay 
JJ  '  77  7J  7y  7j 

where  X  is  a  typical  stress  component  (xx,  yy,  xy)  in  the 
equation  of  motion.  Solving  for  3X/3x  and  3 ^/ay  gives 


a  x  _  i  T a  x  ay  a  x  ay  1 
ax  T  | _TJ  3k  Tic  7j  J 

Ik 

t 

ax  l  f ax  ax  ax  axl 
Ty  "7  [7J  Tic  '  Tic  33  J 


(3) 

(4) 
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where 


T  _  3x  9y  9x  9y 
3j  9T  '  9¥  9j 


(5) 


If  the  k,j  coordinates  assume  discrete  values 
1>  2,  ...  k-1 ,  k,  k+1 ,  ...  k 

max 

1,  2,  ...  j  - 1 ,  j ,  j  + 1 ,  .  . ,  j 


then 


where 


\  =  l^j I lRk|sina  e  = 


max 


*  2Aa  * 


t,  r  ^  J  +  lX  p  ft  9x  9 y 

J  x  9j  ej  \  ~  ex  +  e-j 


Also 


h  *  K 


■  ( 


3_x  9y  9x 


aj  at  '  9¥  If)  ®x  *  ®y  =  J  ^k  *  *y 


(6] 


C7: 


Comparing  Eqs .  (6)  and  (7)  shows  that  a  good  approximation 
to  J  will  be  the  zone  area  (Aa  +  A^)  if 


ex  *  ey 


(8) 


Equation  (8)  is  satisfied  if  the  x,  y  and  k,j  coordinates 
have  the  same  relative  orientation  as  that  shown  in  Fig.  1, 
i.e.,  if  the  unit  vector  obtained 
from  the  ]  x  I  operation  i_.  equal 


to  the  unit  vector  from  e  x  e 

x 
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Fig.  1- -The  x,y  and  k,j 
coordinate  system. 
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Table  1 


Figure  1  also  shows  the  numbering  system  to  be  used 
for  both  nodal  point  and  interior  variables.  Table  1  gives 
the  equivalence  between  the  numbering  system  and  the  k,j 
values . 


Number 

k,j  Interior 

k,j  Exterior 

1 

k-1/2,  j -1/2 

k,j 

2 

k,  j-1 

3 

k-1,  j-1 

4 

k-1,  j 

5 

k-1,  j+1 

6 

k-1/2,  j  +1/2 

k,  j+l 

7 

k+1/2 ,  j  +1/2 

k+l,  j+l 

8 

k+1/2,  j+1/2 

k+1,  j 

9 

k+1,  j-1 

Typical  interior  variables  are  density,  stress,  strain, 
internal  energy  and  area.  Exterior  variables  are  accelera¬ 
tion,  velocity  and  position  vector. 

Equations  (3)  and  (4)  and  Fig.  1  suggest  a  rather 
natural  differencing  scheme;  i.  e., 


1  a_z 

P  dX 


i  y 

7  8  8  1 

P  J  +  P  J 
7  7  8  8 


Wk  + 


z 

6 

,  X,, 

P  J 

+  0  J 

6  6 

1  1 

~T~ 

Z 

y 

8  1  12 

P  J 

+  P  J 

8  8 

1  1 

(1  -  W.  )  - 


E  y  w . 

7  6  6  1  J 

p  iT  +  p  J 

7  7  6  6 


(1 


V 


(9) 
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ft 


I  I 

J 


t 


» 


► 


1  3£_ 
p  3y 


£  x 

7  6  8  1 


£  x 

6  1  14 


pj  +  p  j  k 

7  7  8  8 

PJ 

6  6 

+ 

TT 

i  i 

~T~ 

~T 

£ 

X 

8 

i 

1  2 

P  J 

+ 

P  J 

8  8 

1  1 

£  x  w. 

(1  -  W.  )  +  - a- 6 1- 1 — 3-  - 

v  k"1  p  J  +  p  J 

7  7  6  6 

- 2 - 


(10) 


where  £  =£  -  £  ,  y  =y  -  y  ,  etc.,  and  w,  and  w • 

28788161  K  J 

weight  the  individual  acceleration  components  based  on  the?r 
location  with  respect  to  point  1. 

A  fairly  simple  weighting  scheme  used  in  the  TENSOR 


code 


[8] 


is 


wk = 


8  •  8 

1  4 _ 8_4_ 

8  •  8 

8  4  8  4 


8  •  R 

1  2 _ 6  2 

j  £  .  £ 

6  2  5  2 


W  -  = 


In  Eqs.  (9)  and  (10),  if 


wk  '  wj  1/2 


(11) 


°cJc  *  pdJd 
- 2 - 


pj  +  pJ  +  pj 

1  1  2  2  3  3 

? - 


+  P  J 

4  4 


(12) 


then  the  CRAM  interior  differencing  is  obtained.  Both  HEMP  ^ 
and  CRAM  use  the  same  interior  differencing  scheme.  Equations 
(11)  and  (12)  reduce  the  TENSOR  difference  equations  to  those 
of  CRAM  and  HEMP. 


If  a  k  line  is  to  be  decoupled  from  the  grid,  as 

in  Fig.  2,  then  Eqs.  (9)  and  (10)  permit  this  if  w,  =  1 

+  -  +  ^ 
for  k  and  w^  =  0  for  k  .  For  a  point  on  k  ,  w^  =  1 

and  w^  =  1/2.  Equations  (9),  (10)  and  (12)  may  be  used  to 

write  the  spatial  derivatives  in  Eqs.  (1)  and  (2)  giving 


i 
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-f»  in 

kk* 


w,  =  1 
k 


w.  =0 
k 


Fig.  2--A  decoupled  interior  grid  line.  The  boundary 
stresses  lek*  and  T<y*  do  not  change  when  the  inter 
face  is  viewed  from  below  (w^  =  1)  or  above  (w  =  0)  . 


> 


Fig.  3--The  orthogonal  j",  t.  unit  vectors. 
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1  9xx 

P  dX 


xx  y 

( xx  -  XX*  W  + 

(  xx  -  xx*  \y 

7  8  8  1 

\  7  6/61 

\  8  1/12 

+ 

+ 

L  <f» 

- 

1  3xy 

p  9  y 


™  xy  _x 


7  8  8  1 


I  iXI 

P  9y 


1  9xy 

p  9x 


where 


w  yy.  x 


7  8  8  1 


( xy  -  xyMx  +  ( xy  -  xy*)x 

\  7 _ 6  '  6  1  \  e _  1  /  1 


(yy?  .  ,  (yy<  -  yy*  )x  j 


xy  y 

7  8  8  1 


(xy7  ~  xy:)ysi  +  (xy8  ;  xyt)yi 


(13) 


(14) 


(15) 


(16) 


<J> 


p  J  +  p 

+  _  7  7 


J 

8  8 


(17) 


Along  a  typical  interface  line  joining  two  adjacent 
nodal  points  (Fig.  3)  orthogonal  unit  vectors  5  and  J  are 


k  =  -e  sine  +  e  cos8 

x  y 


j-  =  e  cose  +  e  sine 
x  y 


(18) 


The  stress  components  in  this  coordinate  system  are 

Tele  =  yy  cos2e  +  xx  sin2e  -  2xy  sine  cose  (19) 

FT  =  xy(cos2e  -  sin26)  +  (yy  -  xx)  sine  cose  (20) 

JJ  =  xx  cos2e  +  yy  sin2e  +  2xy  sin6  cose  (21) 

The  acceleration  of  a  point  on  k+  may  be  written  as 

(• 

tt)  =  — r  xx  y  +  xx  y  -  xy  x  -  xy  x  +  kT<*y+ 

Qt  /  (J)+  L  7  8  6  8  2  8  7  8  6  8  2  8  6  6  1 

+  Fk*y+  -  TcJ*x+  -  kT*x+  1  (22) 

y  1  2  J  6  6  1  J 1  1  2  J  v  J 
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» 


» 


f 


Z 


I 


(&)  '  ^  ('yy> 


x0£  '  yy  x  „  +  xy  y  +  xy  y  -  b:*x 

86  828  786  828  661 


rT*x+  -  lc3"*y +  -  Ff*y+ 

1  1  2  6  6  1  1  1  2. 


(23) 


Similarly,  the  acceleration  of  a  point  on  k  may  be  written 
as 


(at)  '  MXV 


+  xx  y  -  xy  x  -  xy  x  -  77*y’ 

64  142  664  142  6  7  6  1 


k¥*y  +  IcJ  x  +  kT*x’  1 
112  661  1  1 2j 


(24) 


~  f-yy  X  -  yy  x  +  xy  y  +  xy  y  +  Et*x" 

<P  L  664  42  664  1  4  2  661 

+  E<*x  +  KT*y  +  Fj"*y~  1 
1  12  *6^61  1  1 2J 


(25) 


These  equations  have  been  derived  assuming  that  Ec*, 

kj*  kk*  and  k j *  are  specified  along  the  6-2  interface, 

6  11  J 

finding  the  corresponding  stresses  in  the  x,y  coordinate 
system,  i. e. , 


xx*  =  7J*  cos2e  +  E<*  sin20  -  2kJ*sin0  cos9  (26) 

yy*  =  IT*  s in2e  +  E<*  cos20  +  2Fpsin0  COS0  (27) 

xy*  "  (IT*  -  Ec*)sin0  COS0  +  7F*(cos20  -  sin20)  (28) 

and  then  substituting  these  expressions  for  xx*,  yy*,  xy* , 
xx*,  yy*,  xy*  into  Eqs .  (13)  through  (16) 

Equations  (22)  through  (25)  will  be  used  to  move  the 
points  on  the  decoupled  grid  line.  They  are  consistent  with 
the  interior  difference  equations. 
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APPENDIX  II 


In  order  to  use  Equations  22  through  25  in  Appendix  I 
the  boundary  stresses  must  be  specified.  In  order  to  solve 
for  these  stresses  we  apply  contact  discontinuity  boundary 
conditions,  which  require  that  the  normal  component  of  stress 

and  normal  component  of  velocity  be  continuous  at  the  bound¬ 
ary  . 

The  unit  vectors  normal  and  tangent  to  the  6-2  inter¬ 
face  may  be  written 

£  =  -sine  ex  +  cose  ey  e^  =  -sine  t+  cosej 

j  =  cose  ex  +  sine  ey  ey  =  cose  t  +  sine j 

The  acceleration  components  on  the  plus  side  of  the  line  may 
be  written 


ak  =  ^rf8k  +  Rc ,  +  R  ^  -  (R  +R  )  TcTc*l 

(J)  L  K  61  7  12  e  f  1  12  J 

ai  =  [g  -  +  R  lg  +  R  Ff  -  (R  +R  )  kj*l 
J  4>  L  J  6  1  7  12b  61  1  2J  J  J 


(2) 

(3) 


The  corresponding  acceleration  components  on  the  minus  side 
are 


a 

k 


(R  +R  )  EE*  -  R  ER- 
6  112  6  16 


(4) 


+ 


+  R  )  EJ*-  R  FT 
12  6  16 


(5) 
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1 


> 


V 


V 


where 


gk  =  '(xx,ay0 ,'xy,  x  )  sine  -  (yy  x  -xy  y  )cose 

K  7881  7881  7881  '  7  8  '  8  1 


(6) 


g,  -  -  (xx  y  -X)  X  )  sine  -  (yy  X  - xy  y  )cos6 

K  61  14  '61  1H;  W/61  14  /61/14; 


st  =  (xx  y  -xy  x  )  cose  -  (yy  x  -xy  y  )si 

J  7881  7881  7881  '  7e'  8T 


ne 


m 


e'j  =  C”t.iy.*'xy,ilc,,)  cose  '  (yy6/lt-xy6iyM)sin° 

y  =  l[p>Vp.Ja]  =  ?kJetp,J,l  (8) 


Equations  2  through  5  have  been  derived  from  Equations 
13  through  16  in  Appendix  I  by  assuming  that  the  boundary 
stress  is  uniform  over  the  entire  6-2  portion  of  the  grid 
line.  These  equations,  therefore,  are  not  completely  con¬ 
sistent  with  the  interior  difference  equations.  The  assump¬ 
tion  concerning  uniform  boundary  stress  is  necessary  since 
the  contact  discont.inuity  boundary  conditions  will  result  in 
an  equation  that  relates  the  normal  components  of  acceleration 
on  the  plus  and  minus  side  of  the  boundary  This  equation 
should  contain  only  one  unknown,  i.e.,  the  normal  component 
of  the  boundary  stress. 

In  order  to  find  the  relation  between  the  above  acceler¬ 
ation  components  at  a  slipping  interface,  we  follow  the  tech¬ 
nique  used  by  Cherry ,  e t  al ,  ^  ^  and  attach  a  coordinate  system 
to  the  point  on  the  minus  side,  with  t  and  J  being  the  unit 
vectors  normal  and  tangent  to  the  slip  line  at  the  point  to 
be  moved. 
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r 


? 


7 


11  v  and  v  are  the  velocities  on  the  plus  and  minus 
side  of  the  slip  line  then,  from  Equation  1,  the  normal  com¬ 
ponents  of  velocity  are  given  by 


V 

II 

• 

•  + 
-x 

sin0 

•  + 

+  y 

COS0 

V 

11 

• 

•  — 

-x 

sin0 

*  y- 

COS0 

If 


,  -  _  3x_L  -  9y  ~ 

3£  y£  9£ 


then 


R, 


1/2 


sin0 


cos8 


CIO) 


Substituting  Equation  10  into  9,  and  equating  normal  velocity 
components  gives 


R£v  *k  =  R£v  •  k 


Cll) 


or 


7 


f 


-xyl+y‘*l  -  -*+y~!L+y+xl  (12) 

Since  the  1<,£  coordinate  system  is  attached  to  the  minus  side 
of  the  slip  line  then 


<L  , 

dt  '■ 


x  y£+y 


1  9  r  ’ +  ~  +  ~  \ 

)  =  3t  C-x  yt*y 


k,£  constant 


(13) 
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Tlie  left  side  of  Equation  13  may  be  written 

-n 

(14) 

=  R^ak  +  (v  4)Cv;.J)  -  (v-.t)(^4) 

Since 


d 

3t 


(VM,  .  f^vM) 


+ 


9 

•3k 


and 

.  + 

iv  -v) *k  =  0 

Then  the  right  side  of  Equation  13  may  be  written 

(15) 

(v  .*)  (v^£)  -  (v+-v>£v+4 


O  f  •  *  -  •  +  - 

3t('x  y%*y  xO 

*  +  (  v“  •  j)  - 


In  Equations  14  and  15  a^  and  a^  are  given  by 

a+  =  -  dl*+)  Z±  +  iC/.!)  VJl 

k  3t  R^  dt  R^ 


d(x~ )  y l  +  d (y  )  ^ l 
dt  R£  dt  R^ 
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and  arc  the  sane  acceleration  components  given  by  Equations  2 
and  4.  We  have  also  used  the  relation 


-X  y£+y  x£  =  (v“ • k)  (v£  •  j )  -  (v"*X)(v£-k) 


(16) 


in  order  to  derive  Equations  14  and  15. 

Substituting  Equations  14  and  15  into  Equation  13 

gives 


Vak  =  Ac 


(16) 


where 


A  = 
c 


(v  -V  )  *X(v£+v£)  •  k 


R, 


(17) 


Solving  Equation  16  for  IcF*  gives 


=  ^k^X  +  -  »VA( 

(18) 


(R  +R  )(4>  +4") 
6  1  12 


Equation  18  relates  the  normal  component  of  stress  (Iclc*) 
at  the  slipping  interface  to  interior  zone  variables.  In  order 

to  move  the  boundary  points  then  kj*  in  Equations  22  through  24 
in  Appendix  I  must  also  be  specified. 

For  a  tied  point  Aq  in  Equation  18  is  set  equal  to  zero, 

since 


(v  -X  ) • £  =  0 


Also  the  tangential  stress  component,  for  a  tied  point,  is 
obtained  from  Equations  3  and  5.  Since 


81 


I 


R- 1759 


aj  =  ai 


then  for  a  tied  point 


<f>'g*-4>+g  •  +  4-  (R  FT  +  R  TcJ*  )  +  <{>"  (R  Fj"  +R  FT  ) 

^  &]  bJ  61  j7  12  Je  61  J6  12 


Ff*  = 


(R  +  R  )(<p 

6  1  12 


v 


^  ► 
4  t 


(19) 
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Plastic  flow  is  due  to  the  inability  of  real  materials 
to  support  unlimited  values  of  shear  stress.  In  the  code  the 
deviatoric  stress  components  are  modified  such  that  the  re¬ 
sulting  stress  state  is  consistent  with  a  Mises  yield 
criterion. 

If  the  second  deviatoric  invariant  (J)  is  greater  than 
a  specified  value  (1/3  Y2),  then 


CD 


where  is  the  adjusted  stress  deviator 

A 

is  the  stress  deviator  calculated  by  assuming 
that  the  total  strain  rate  is  elastic,  and 


(?) 


For  a  triaxial  test,  Y  corresponds  to  the  maximum  allowable 
stress  difference  at  failure. 

Rupture  initiation  is  being  modeled  by  accumulating 
the  difference  between  /TJ  and  Y  during  yielding.  When 
this  accumulation  reaches  a  specified  value  then  the  point 
at  the  fault  surface  enters  the  slip  routine.  Between  two 
consecutive  cycles,  n  and  n+1,  the  accumulation  takes  the 
form 


83 


» 


R- 1 759 


» 


y 


4 


? 


V 


T 


o 


3 


Rupture  occurs  if 


n+1 

e 


>  W 


(4) 


where  W  is  a  specified  function  of  distance  from  the 
initial  point  of  rupture  (the  focus). 


Equation  (3)  is  similar  to 
where  the  plastic  work  (En+^)  is 


r.n+1  T,n  ,  Y2  /IT  -  Y 
E  =  E  +  ^  — ? — 


a  plastic  work  criterion, 
given  by 


(5) 


Equations  (3)  and  (5)  differ  only  by  the  factor  Y2/3y. 

We  have  been  successful  in  both  controlling  rupture 
velocity  and  reducing  the  stress  concentrations  at  the  end 
of  the  fault  by  allowing  W,  in  Eq.  (4),  to  be  a  specified 
function  of  distance  from  the  point  of  rupture  initiation. 
The  functional  form  that  has  been  used  is 


W 


6c 

c 

7 


x  +  L/2\2["l  1  x  +  L/2] 

; — I?  '  7  — a — _ 

+  3  £-y^] 


0  <  x  + 


(6a) 

(6b) 


where  L,  c  and  d  are  input  parameters.  The  rupture  is 
constrained  to  lie  between  -  L/2  <_  x  <_  L/2  where  L  is 
the  fault  length. 
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